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AN  ALTERNATE  APPROACH  TO  OPTIMUM 
BEARING  ESTIMATION 

INTRODUCTION 


rhe  optimum  passive  bearing  estimation  problem  is  examined  in  this 
report  from  a 'k-w  space"  point  of  view.  (The  two-dimensional  Fourier  trans- 
orm  of  the  acoustic  field  in  the  space-time  aperture  generates  the  energy 
distribution  in  the  k-w  domain.)  Since  the  representation  of  the  acoustic  field 
m k-  w space  has  proved  useful  in  a number  of  applications  and  since  recent 
sonar  systems  essentially  generate  a sampled  function  of  the  energy  distribu- 
tion in  k-w  space,  this  viewpoint  is  expected  to  provide  additional  insight  into 
the  processing  required  for  optimum  bearing  estimation. 

particular,  guided  by  an  intermediate  result  of  MacDonald  and 
Schultheiss  and  a result  by  Bouvier  and  Ianniello,  2 it  is  postulated  that  the 
maximum  likelihood  estimator  (MLE)  merely  determines  an  estimate  of  the 
location  of  the  primary  peaks  of  the  energy  distribution  in  k-  oj  space  and  com- 
buies  these  estimates,  through  a linear  combination,  to  arrive  at  an  overall 

eSt“nate  of  the  target  bearing.  This  postulate  is  examined  in  this  report 
and  the  performance  of  such  an  estimator  is  compared  with  the  Cramer-Rao  lower 
bound  under  the  same  conditions  assumed  by  MacDonald  and  Schultheiss. 

PRELIMINARIES 


„ ..  ConsicJer  a finite,  one-dimensional  array  to  be  described  by  an 
effective  sensor  density  p(x),  where  p(x)  includes  the  effect  of  sensor  location 
and  sensor  shading  and  is  considered  to  be  zero  for  values  of  lx  I > L/2  where 
L is  the  total  length  of  the  array.  In  practice,  p(x)  is  usually  discrete,  ’but 
mathematically  it  may  be  treated  as  a continuous  function  of  x multiplied  by  an 
appropriate  sampling  function.  If  h(x,  t)  represents  the  acoustic  field  at  any 
point  x and  time  t,  then  the  output  from  a conventional  time  delay  beamformer 
steered  to  an  angle  8 with  respect  to  broadside  is 


b ( 4 , t)  = / dx  h(x,  t + £ x)  p(x) , 

• — CO 

where  | - - and  c is  the  speed  of  sound  in  the  medium. 


(1) 


over  a 
output 


Furthermore,  consider  that  this  beam  output  is  spectrum  analyzed 
sliding  temporal  window  of  duration  T.  The  spectrum  analyzed  beam 
is  then  described  by 


1 


TR  5439 


B = 1/T  / ^ dtb  (|  , t)  e~i2n  ft.  (2) 

t'  - T/2 

Substituting  (1)  into  (2)  and  assuming*  that  T>  > L/c  yields 

B = 1/T  //  dx  dt  h(x,  t)  p(x)  II  (^f~j  e _l27rft  e i2?rf|x  ^ 

where  II  ^ is  the  familiar  rectangle  function  centered  about  t'  and  of  unit 
height  and  duration  T. 

For  comparison  purposes,  the  direct  two-dimensional  Fourier  trans- 
form of  the  acoustic  field  h(x,t),  as  seen  through  a spatial  aperture  described 

r 

by  p(x)  and  a temporal  aperture  described  by  1/T  is 

8=  1/T  //  dx  dt  h(x,  t)  p(x)  II  c-12^  e _i2tt^X>(4) 

— CO  1 

which  represents  the  acoustic  field  in  k-oj  space,  that  is,  at  a temporal 

frequency  w = 2 n f and  a spatial  frequency  k =£«  = ^ SI^  . Comparison  of 
(3)  and  (4)  shows  that  B(k,  w ) = 8 (-k,  w ).  Thus,  spectrum  analyzing  the 
various  beams  will  generate  the  acoustic  field  in  k-w  space,  which  merely 
reaffirms  the  fact  that  time  delay  beamforming  is  equivalent  to  performing  a 
spatial  Fourier  transform. 

In  practice,  (3)  must  be  modified  by  an  appropriate  sampling  function 
since  beams  are  formed  only  at  discrete  values  of  t,  and  are  spectrally 
analyzed  only  at  discrete  temporal  frequencies.  For  the  initial  analysis, 
however,  the  spectrum  analyzed  beam  output  B ( £ , f)  is  considered  to  be  a 
continuous  function  of  | and  f.  Appropriate  modifications  are  made  later. 

In  general,  the  acoustic  field  h(x,  t)  will  consist  of  a signal  and 
noise  component.  If  the  signal  is  a plane  wave  arriving  at  an  angle  0O with 
respect  to  the  normal  of  the  x-axis,  as  shown  in  figure  1,  h(x,t)  may  be 
expressed  as 

h (x,  t)  = s (t-£0x)  + n (x,  t), 

where  £0=  s^6- . Thus,  provided  the  noise  exhibits  no  spatial  structure,  the 

♦This  assumption  implies  that  the  observation  time  T is  much  greater  than 
the  travel  time  of  the  signal  wavefront  across  the  array. 
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acoustic  energy  distribution  y(  £ , f)  = | B(  £ , f)  | 2 will  exhibit  a distinct 
maximum  when£  = £0  . In  accordance  with  the  postulate  mentioned  earlier, 
it  is  desired  to  obtain  the  mean  and  variance  of  the  estimate  | 0 at  which 
y is  a maximum. 


ARRAY  NORMAL 


INCIDENT  SIGNAL  WAVE 


/ SENSOR  LOCATION 

Figure  1.  Array  Geometry 


MEAN  AND  VARIANCE  OF  ESTIMATING  THE  LOCATION  OF  THE 
GLOBAL  PEAK  OF  ENERGY  DISTRIBUTION  IN  k-w  SPACE 

Consider  the  spectrum  analyzed  beam  output  B ( £ , f),as  given  by 
(3),  to  be  a continuous  function  of  £ and  f.  As  shown  in  figure  2,  the  proposed 
processor  searches  for  the  global  peak  of  the  square  of  the  magnitude  of  B. 


B(fcf) 

2 

y<e,f) 

SEARCH 

rUR 

PEAK 

Figure  1.  Proposed  Processing  Scheme 

A 

Accordingly,  at  a given  temporal  frequency  f,  the  value  of  £ =£0  at  which  y is 
a peak  is  considered  to  be  an  estimate  of  the  target  location.  The  mean  and 
variance  of  this  estimate  are  given  by 


varUo>  =i 


BE  (y') 

*4 


E u ) = t - E (y’> 

Uo)  \)  3E  (yf) 


* = 40 


* = «o 


where  y*  = — y ( £ , f).  Since  taking  the  expectation  is  a linear  operation,  it 
d4  3E  (v') 

may  readily  be  shown  that = E (y  "),  so  that  the  above  relations 

become 
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1 


A 

var  (So) 


E (y'2)  - E2  (y») 
E2“<y"> 


S=S( 


E (So)  = «o 


E (y") 


.(6) 

S = So 


The  various  expectations  required  above  shall  be  computed  under  the  following 
assumptions: 

1.  Signal  and  noise  are  spatially  homogeneous  and  temporally 
stationary,  zero  mean,  independent  Gaussian  processes. 

2.  The  observation  time  T is  large  with  respect  to  the  travel  time  of 
the  signal  wavefront  across  the  array. 

3.  The  signal  wavefront  may  be  regarded  as  a plane  wave. 

Later  assumptions  shall  include: 

4.  The  noise  is  spatially  incoherent  and  uniform. 

5.  The  effective  sensor  density  is  discrete  and  uniform  (no  shading). 

6.  The  observation  time  is  large  with  respect  to  the  signal  and  noise 
correlation  times. 


Throughout  this  analysis,  one  will  require  the  crosscorrelation 
between  two  separate  points  in  the  space-time  aperture.  The  assumption  of 
spatial  homogeneity  and  temporal  stationarity  allows  this  correlation  to  be 
expressed  as  a function  of  spatial  and  temporal  differences;  that  is, 

E [h  (Xj,  tj)  h (x2,  t2)J  = Rh  (x,r  ) , (7) 

where  x = x-  - x and  r = t - t . 

Assumption  3 allows  the  further  simplification: 

Rh  (x, r ) = Rg  (r-£ox)  + \ (x,  t),  (8) 

where  R and  R are  the  signal  and  noise  correlation  functions,  respectively, 
s n 

Notice  that  Rg  depends  only  on  the  difference  (t  -£0x).  Finally,  assumption 
4 (uncorrelated  noise)  allows  R^  (x,  t ) to  be  expressed  as 

\ (x,  r ) = Rg  (t-|qx)  + 6 (x)  Rn  ( r ),  (9) 


1 

} 

* 4 

*■  ‘i 

) 
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where  8{x)  is  the  impulse  function  and  R (t)  refers  to  the  noise  autocorrela- 
11 r- 11 i n 


tion  function  at  each  sensor. 


The  derivation  of  the  expectations  required  in  equation  (6),  which  is 
relatively  straightforward  but  tedious,  is  presented  in  appendix  A.  The  results 
are  quite  general— being  subject  only  to  the  first  two  assumptions  listed  above. 
In  particular,  it  is  shown  that,  provided  f > 1/T  the  expectation  and  variance 
of|o  is 


E (?o>  = €o 


//  dxdT0(p  * p)  (icux) 


//  dxdr0(p  * p)  (iojx)2 
w2  [//dxdT0(p  * p)  x2]2 


(10) 


4=£o 


■{  2 Re 


|//dxd  t 0 (p*xp)|  2 - 2 (j[ dxd T0(p  * xp) 


(11) 


'ff  dxdT0(p  * p)  ff  dxdT0(xp  * xp)  - ^//dxdT0x  (p  * p)^  2 } 


where 


0 = 0 (x,  t ) = Rh(x,  r)  A ( t/T)  e-1UJTe1C0  ^ X 

A(t/T)  = ( 1-lTl 
\ 0 


(12) 


(Triangle  Fimction) 


1-lTl/T  I rl<  T 

otherwise 
Re(  ) denotes  the  real  part  of  ( ) 

* denotes  correlation;  e.  g. , g (x)  * u(x)  = / g(t)  u(x  +- 1)  dt 
T = observation  time 


R^  (x, r ) = E[h(x1,  tj)  h(x1  - x,  t -T)]  = crosscorrelation 


between  the  acoustic  field  at  two  separate  points  in  the  space- 
time  aperture. 


The  integrals  are  to  be  evaluated  at  the  value  of  £ corresponding  to 
the  true  target  bearing  0Q  and,  unless  otherwise  noted,  all  integrations  arc 
from  - co  to  + ® . 


These  expressions  show  the  effect  of  an  arbitrary  spatial  weighting 
function  p(x)  on  estimating  the  location  of  the  peak  of  the  energy  distribution 
at  a given  temporal  frequency  f > 1/T.  As  indicated  in  appendix  A,  the 
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expressions  may  also  be  modified  to  include  any  arbitrary  temporal  weighting 
function  w ( t ) by  simply  replacing  A ( t/T)  in  (12)  by  1/T  [w(  r ) * w(  r )] 
and  restricting  temporal  frequencies  to  those  greater  than  the  width  of  the 
corresponding  spectral  window. 


To  gain  further  insight  into  the  performance  of  this  estimator,  it 
shall  now  be  assumed  that  the  signal  is  a plane  wave  and  that  the  noise  is 
spatially  uncorrelated  and  uniform.  As  shown  in  appendix  B,  these  assump- 
tions allow  each  double  integral  in  (10)  and  (11)  to  be  expressed  as  a sum  of  the 
product  of  two  single  integrals.  Furthermore,  the  estimator  becomes  unbiased 
since,  under  these  conditions, 

// dxdr0(p  * p)  x =0. 

The  result,  as  obtained  from  appendix  B,  is  that  E (?)  = ^Qand 


var  (?)  = — — 

w2  [/  x2  (p  * p)  dx]“ 


• { N/S[2 /pdx  f x2  (p*  p2)  dx  - /p2dx /x2  (p*  p)  dx] 
+(n/S^2  / x2  (p2*  p")  dx  }, 


where 


S - S(f)  - / dQ-Gs(o)  sine2  T (f-o)  = Gg  (f)%  sine2  Tf 

N = N(f)  = / da'Gn  (o)  sine2  T (f-o)  = Gr  (f)*  sine2  Tf 
Gg  (f)  is  the  signal  spectral  density 
Gn  (f)  is  the  noise  spectral  density. 

Note  that  one  may  regard  S(f)  and  N(f)  as  being  the  signal  and  noise 
power,  respectively,  within  the  analysis  band  which  is  centered  on  frequency  f. 


Thus  far,  no  particular  restriction  has  been  placed  on  the  sensor 
density  p(x).  However,  the  uncorrelated  noise  condition  (which  has  been 
applied)  is  rather  unrealistic  unless  p(x)  is  discrete.  If  it  is  further  assumed 
that  p(x)  is  uniform  (no  shading),  then  p2  in  (13)  may  be  replaced  by  p.  The 
variance  then  reduces  to 


<!o>  = 


2 r 2 

w L/x  (p*  p)  dxj 


[(N/S)/pdx  + (N/S)  2j 
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Mathematically,  a uniform,  discrete  array  consisting  of  M sensors  may  be 
described  by  a density 


where  6(x)  is  the  impulse  function  and  x.  is  the  coordinate  of  the  jth  sensor 


Substituting  (15)  into  (14)  yields 


This  expression  represents  the  uncertainty,  under  the  first  five 
assumptions  listed  previously,  of  locating  the  global  peak  of  the  energy  dis 
tribution  in  k-o;  space  at  a single  temporal  frequency  u>.  Since  = sin  0O, 

c 

it  may  readily  be  shown  that 


where  6Q  is  the  estimate  of  the  signal  bearing  angle  with  respect  to  broadside 
Thus,  the  variance  of  e is 


In  these  expressions,  notice  that  the  effect  of  the  signal  and  noise 
spectrum  is  included  through  the  relationships 


S (f)  = G (f)  * sine  Tf 


N (f)  = G (f)  * sine  Tf 


In  particular,  if  the  observation  time  is  large  with  respect  to  the  signal  and 
noise  correlation  times  (assumption  (>),  that  is,  if  G and  G are  broad  with 
2 s n 

respect  to  sine  Tf,  then  S (f)  »1/T  G (f)  and  N(f)^  1/T  G (f).  Substituting 
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these  spectral  densities  into  (18)  and  comparing  the  result  with  expression  (16) 

of  MacDonald  and  Schultheiss1  shows  that  the  two  expressions  are  equivalent 
at  a single  temporal  frequency.  Therefore,  at  a single  temporal  frequency, 
the  procedure  of  forming  an  infinite  number  of  beams  and  searching  for  the 
peak  of  the  energy  distribution  is  optimum  in  the  sense  that  the  uncertainty  of 
this  estimate  is  equal  to  the  Cramer-Rao  lower  bound  (CRLB). 


The  question  of  combining  the  bearing  estimates  obtained  at  other 
temporal  frequencies  shall  now  be  discussed.  For  this  purpose,  assume  that 
each  of  the  beams  is  spectrum  analyzed  according  to  (2)  at  evenly  spaced, 
discrete  frequencies  f = n/T  with  n = 1,  2,  3,  . . . , N^.  The  reason  for 

this  choice  is  that  if  the  temporal  periodicity  of  the  acoustic  field  is  matched 
to  the  observation  time  T,  then  the  minimum  sampling  interval  required  to 
properly  describe  the  acoustic  field  is  l/T. 


Consider  now  that  the  peak  of  the  energy  distribution  is  estimated  at 
each  temporal  frequency  f . According  to  the  previous  analysis,  each  estimate 

n g 0 

| may  be  considered  a random  variable  whose  mean  is  and 

2 

whose  variance  ° is  given  by  (16)  with  f replaced  by  f . Furthermore,  if  T 
n n 

is  large  with  respect  to  the  signal  and  noise  correlation  times,  it  is  reasonable 
to  assume  that  each  estimate  (£  n)  is  statistically  independent. 


A natural  scheme  for  combining  these  estimates  is  to  form  a new 


random  variable  £ 
let 


eq. 


which  consists  of  a linear  combination  of  % ; that  is 

- n 

N 

o 

= I * ^ 


eq. 


n=l 


In 


(20) 


where  N represents  the  number  of  independent  estimates  available  and  a is 

a weighting  factor.  Under  the  constraint  that  the  mean  of  £ also  be  equal 

to  £,  and  that  its  variance  be  a minimum , it  may  be  shown  that  £ should  be 
bo  eq. 

formed  according  to 

N 

f */%  *n,  (21) 

n = 1 
n=l 

and  that  the  variance  of  £ is 

^eq- 


^eq.  No 


l/c 


8 


Relation  (21)  shows  that  each  estimate |n  should  be  weighted  by  the  inverse  of 
2 

its  variance  o . Referring  to  (16),  this  implies  that  estimates  having  a 
large  ratio, 


should  be  weighted  more  heavily 


Substituting  (16)  into  (22)  and  using  (17)  and  the  fact  that  S (f^) 
G (f  ) and  N(f  )*1/TG  (f  ) one  obtains 


where  0 is  the  best  estimate  of  the  target  bearing 


= the  highest  frequency  processed, 
o 

The  above  expression  is  identical  to  the  CRLB  derived  by  MacDonald 
and  Schultheissl  (their  equation  16)  and  agrees  with  Carter'*  when  M=2.  It 
follows  that,  under  the  sLx  conditions  listed  earlier,  the  CRLB  can  be  achieved 
by  a processor  that  estimates  the  angle  6 of  the  peak  of  the  energy  distribu- 
tion in  k-w  space  at  each  processed  frequency  f = n/T,  n = 1,  2,  . . . , N 

and  obtains  a weighted  sum  of  the  N(  bearing  estimates.  The  particular 
weighting  to  be  used  in  the  averaging  is  proportional  to  the  inverse  of  the 
uncertainty  associated  with  each  estimate. 

Thus  far  it  has  been  assumed  that  a continuous  beamformer  is  avail- 
able for  estimating  £ . Obviously,  in  practice,  only  a finite  discrete  number 
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of  beams  can  be  formed.  This  dilemma  leads  directly  into  the  classical  inter- 
polation problem.  One  practical  solution  is  to  form  a "large  number"  of  beams 
and  to  use  a simple  interpolation  scheme,  such  as  3-point  interpolation,  in 
order  to  estimate  the  location  of  the  global  peak.  Since  the  beam  pattern  is 
known,  another  approach  is  to  crosscorrelate  the  finite  number  of  spectrum 
analyzed  beam  outputs  with  the  amplitude  pattern  of  the  sensor  density.  This 
crosscorrelation  must  be  performed  at  each  processed  temporal  frequency 
and  must  precede  the  squaring  operation  shown  in  figure  2.  Obviously,  the 
amount  of  processing  required  is  quite  extensive,  but  it  is  possible  to  com- 
pletely recover  the  continuous  beam  output  if  the  original,  finite  beams  are 
formed  at  the  Nyquist  rate. 

The  interpolation  problem  serves  to  highlight  the  fact  that  the  CRLB 
can  only  be  approached  by  any  practical  estimator. 


APPLICATION  OF  BEARING  ESTIMATION  SCHEME  TO 
RECENT  SONAR  SYSTEMS 

The  "k-u;  space"  concept  has  led  to  an  alternate  processing  scheme 
which,  under  the  condition  of  a broad  signal  spectrum,  is  capable  of  achieving 
the  CRLB.  In  view  of  the  large  amount  of  required  processing  and  the  fact 
that  existing  split -beam  correlators  are  nearly  optimum,  this  processing 
scheme  is  not  recommended  for  the  sole  purpose  of  estimating  the  bearing  of 
a broadband  signal.  However,  the  technique  is  quite  attractive  for  recent 
sonar  systems,  particularly  since  the  input  quantities  required  for  this 
estimator,  e.  g.,  spectrum  analyzed  beam  outputs,  are  already  available. 

When  the  signal  spectrum  becomes  more  impulsive  (that  is,  when 
one  can  no  longer  assume  that  the  observation  time  T is  large  with  respect 
to  the  simal  autocorrelation  time),  (16)  and  (19)  show  (with  f replaced  by 
fR)  ^at  “ie  variance  of  | becomes  unduly  large,  because  of  the  reduced 

signal-to-noise  ratio,  unless  the  processing  frequency  f is  aligned  on  the 

n 

particular  signal  tonal  frequency  f^.  Hence,  additional  frequency  tracking  is 

required.  This  leads  directly  to  the  concept  that  bearing  estimation  is 
actually  a two-dimensional  search  (in  k and  w ) for  the  primary  peaks  of  the 
energy  distribution.  If  the  temporal  frequency  rate  is  sufficiently  high,  this 
frequency  tracking  may  take  the  form  of  simple  interpolation  between  fre- 
quency bins.  Notice  that  the  variance  of  the  resultant  bearing  estimate  e at 

the  nominal  frequency  fn  should  correspond  to  that  obtained  by  DIFAR*-like 
automatic  target  followers  (AT  F)  operating  on  a single  signal  tonal  frequency. 

As  before,  if  the  frequency  fn  associated  with  each  estimate  is  sufficiently 
separated  so  that  the  £ n are  independent,  the  estimates  may  be  combined 
"across  frequency,"  in  the  manner  indicated  by  (21),  to  obtain  an 

‘Direction  Finding  And  Hanging. 
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overall  ,fbest"  estimate  whose  variance  is  given  by  (22),  where  now  Nq 
represents  the  number  of  independent  estimates. 

From  the  previous  comments,  it  is  clear  that  the  estimator  is  also 
attractive  from  the  viewpoint  of  replacing  the  conventional  DIFAR-like  ATF 
because 

1.  The  estimator  operates  in  an  open  loop  and  is  not  afflicted  by 
the  acquisition  and  gate  widths  associated  with  closed  loop 
trackers. 

2.  The  estimation  scheme  combines  bearing  estimates  obtained 
from  multiple  tonal  frequencies.  The  estimator  is  therefore 
more  immune  to  fading  tonals  than  is  the  conventional  ATF. 

3.  The  tracking  threshold  equals  the  detection  threshold  since  the 
estimation  scheme  merely  locates  the  peak  (in  frequency  and 
bearing)  of  detected  tonals. 

The  effects  of  multiple  targets  and  frequency  and  bearing  dynamics 
have, thus  far,  been  purposely  avoided.  It  is  clear,  however,  that  although 
the  basic  bearing  estimation  procedure  remains  the  same^  the  effects  give 
rise  to  a sorting  problem.  Before  the  various  estimates  ( or  0n)  can  be 

combined  "across  frequency,  " one  must  ensure  that  they  are,  in  fact,  from 
the  same  target. 


CONCLUSIONS 

It  has  been  shown  that  the  processing  scheme  suitable  for  estimating 
the  bearing  of  a plane  wave  in  spatially  incoherent  noise  is  one  that  locates  and 
linearly  combines  the  estimates  of  each  global  peak  of  the  energy  distribution 
in  k-ui  space.  For  sonar  systems  that  spectrum  analyze  each  of  a finite  num- 
ber of  preformed  beams,  this  amounts  to  estimating,  at  each  processed 
frequency,  the  value  of  sin  B at  wdiich  the  magnitude  squared  of  the  beam 
outputs  is  a maximum.  These  estimates  may  then  be  combined  by  forming  a 
weighted  average  where  the  weighting  is  proportional  to  the  inverse  of  the 
uncertainty  associated  wdth  each  estimate.  It  has  been  showai  that,  for  broad 
signal  spectra,  such  an  estimator  is  capable  of  achieving  the  Cramer-Rao 
lower  bound  and  is,  therefore,  optimum. 

For  impulse  signal  spectra,  although  a similar  procedure  is  indicated, 
an  analytical  expression  for  the  variance  of  the  overall  estimate  is  difficult 
to  obtain  — partly  because  of  the  unknown  statistics  of  adjacent  spectral 
estimates. 
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Even  f°r  general  use  with  multiple  targets  and  frequency  dynamics 
the  bearing  estimation  scheme  is  considerably  complicated  by  the  necessary 
sorting  problem.  J 
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DERIVATION  OF  THE  EXPECTATIONS  AND  VARIANCES 


In  this  appendix  we  wish  to  derive  the  expectation  and  variance  of  the 
estimated  coordinate  of  the  global  peak  in  y = I Ii| 2 , where  B is  the  spectrum 
analyzed  beam  output  as  expressed  by  (3)  in  the  text,  which  is 


For  this  purpose,  only  the  following  two  assumptions  are  made: 

1.  Signal  and  noise  are  spatially  homogeneous  and  temporally  stationary 
zero  mean,  Gaussian  processes. 


The  observation  time  T is  large  with  respect  to  the  travel  time  of  the 
signal  wavefront  across  the  array.  (This  assumption  has  already  been 
made  in  expression  (A-l). ) 


In  accordance  with  (6)  of  the  text,  we  need  to  evaluate  the  expectation 
of  derivatives  of  y as  well  as  the  expectation  of  (y')^.  The  expectation  of 
derivatives  of  y may  be  determined  in  the  following  direct  manner.  From  (A-l), 

y = BB*  = l/T2  ////  dx1dx2  d^dt,,  hfXj.tj)  h(x2,t2)  p(x1)  p(x2) 


where  cj  = 2 rr  f and  all  integrations  are  from  - =o  to  + 11 . Now  the  nth  deriva 
tive  of  y with  respect  to  £ is 


and  the  expectation  of  the  nth  derivative  of  y is 
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E [y(n)]  = 1/T2  ////dx1dx2dt1dt2  E [h(xlf  tj)  h(x2,t2)]p(Xl)  p(x2) 

(tj-t*)  (t2-f)  „ iu4(Xl-x2) 

• n n [iw(x..-x  )]  e e • (A-4) 

T T 1 ^ 

Because  of  the  spatial  homogeneity  and  temporal  stationarity  of  the  acoustic 
field  h(x,t),  we  may  write 

E [h(x1>t1)  h(x2,t2)  ] = Rh  (x1-x2,  tt-t2)  , (A-5) 

and  introduce  the  following  integration  variables  in  (A-4): 

x = x -x  in  the  integration  over  x 
12  ^ 

t = t -t  in  the  integration  over  t . 

12  c 


Then,  (A-4)  becomes 

r (n)-.  2 

E [y  ]=  1/T  //// dx,  dxdt,  dTRh  (x,  r ) pfx^  p^-x) 


(t,  -f)  (tj-T-f) 


„ * _ 1 >n  -i«T  icufx 

• n II  (iujx)  e e 


T T 

If  we  perform  the  integration  over  x^  and  and  note  that 

/ dx1  p (x^  p(Xl-x)  = p(x)  * p(x) 


and 

(t^f)  (tj-t'-r) 

/dt,  n II =n  (T/T)  *H(T/T)  = T a(t/T), 

1 T T 

I 1-  M/T  MiT 

where  A(  t/T)  = \ 

( 0 otherwise 


(A-6) 


(A-7) 


(A-8) 


and  that  * denotes  correlation,  (A-(i)  reduces  to 
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E [y^]  ~ 1/T  //dxdr  6 (x,  r)  (iwx)n[p(x)  * p(x)]  , (A-9) 

where 


4><x,r ) = Rh(x,T  ) a(t/T)  e~luJT  e1"4  X . 


(A-10) 


The  evaluation  of  E(y'2)  is  more  readily  accomplished  in  an  indirect 
manner.  Thus. 


= I bI  2 = 


BB* 


and 


<y 

where  Re  ( 
becomes 


•>2^ 


= (BB*'  + B*B')2  = 
) denotes  the  real  part  of  ( 


2 Re(B*B'B*B'  + BB*B'B'*)  , 

).  The  expectation  of  (y')2  then 


E(y'2)  = 2 Re  [E(B*B'B*B')  + E(BB*  B’B'*)]  . (A 

Since  the  spectrum  analyzed  beam  output  B is  also  a zero  mean  Gaussian 
process,  we  may  use  the  familiar  "chain  rule”  for  the  expectation  of  four 
products  to  obtain 


E(y'2)  = 2 Re{2  [E(B*B')]“  + [e  (BB)]  E(B'B’) 


+ E(BB*)  E(B'B'*)  + 


E (BB') 


E(B*B')2|  | 


(A-12) 


In  a manner  completely  analogous  to  the  derivation  of  E fy^] 
it  may  be  shown  that 

E [B^B^]  = 1/T  //dxd  r<t)(x,r)[(-iwx)n  p(x)  * (iwx)kp(x)]  (A-13) 
and 


-i2wt' 

E[B^k)B(n)]=  — //  dxd  t ^(x,  T)[(iu.x)n  p(x)  elw^X  (iwx)k  p(x)ei>' - X], 

(A-14) 

where 

v(x,  r)“Rh(x,  r)  A(  r/T)  sine  2f  (T-|t|)  (A-15) 
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sine  Q = 


B^  = kth  derivative  of  B with  respect  to  | . 

Expressions  (A-13)  and  (A-14)  allow  the  evaluation  of  each  of  the 
expectations  in  (A-12).  For  example,  E(B'B'*)  is  obtained  from  (A-13)  with 
n=k=l  and  E(BB')  is  obtained  from  (A-14)  with  k=0,  n=l.  The  complete  ex- 
pression for  E(y'2),  in  shorthand  notation,  is 


E(y’2)  = — “ Re  { | //  dxdT(j>(p  * xp)| 

T2 

- | //  dxd  r (peiw-x  * x p e1  w ^ x)  | 2 

+ SI  dxdT<j>(p  * p)//  dxdT4>(xp  * xp)  (A-16) 

rr  , , , -iw£x  -iw£x,  /•/«,  , . iuj£x  iwtx 

- JJ  dxdry(pe  s(c  pe  )// dxdnji(xpe  sf:  xpe  ) 

- 2[//  dxdr<)>(p  * xp)  ] 2 } . 

This  relationship  may  be  simplified  considerably  via  the  following 
argument.  The  second  term  of  (A-16)  may  be  expressed  as 

I fr,  , iwfcx  . i2u>£x  . I 2 
I JJ  dxdT  ^e  * (pe  * xp)| 


Then  using  (A -15)  and  (A-10),  this  may  be  written  as 


i2  u>|  x 


I r • (pe  * * xp)"]) 

|//dxdT  (b(p  xp)  I elu>T  sine  2f  (T- Irl ) J|  • 

L (p*xp) 

Now  the  magnitude  of  the  term  in  the  brackets  [ ] is  always  much  less  than 
one  provided  f > 1/T.  Thus,  provided  we  use  temporal  frequencies  in 
excess  of  l/T, 

|//dxdr<t>(p*  xp)  2»  |j/dxdr^(peiw^X  * xpe1"'*)  J 2 . 

Similarly,  the  third  term  of  (A-16)  will  always  be  much  greater  than  the 
fourth  term  if  f > 1/T.  Hence,  under  this  minor  restriction  of  temporal 
frequencies,  (A-16)  reduces  to 


E (y’2)  = 


(A-17) 


~Y  Re 1 1 Jfdxd r<j>(p  * xp)  |2  - 2 [//dxdr4>p  * xp)]2 
+ //  dxdr  (P  * P)  * JJ  dxd  r <i>  (xp  * xp)  j> , 

and  from  (A-9), 

E(y')  = ^r//dxdr4>x(p  * p)  (A-18) 


and 

_ 2 ,,  2 

E(y")  = — — JJdxdr  <b  x (p  * p).  (A-19) 

T 

Finally,  the  substitution  of  (A-17)  through  (A-19)  into  (6)  of  the  text 
leads  to  expressions  (10)  and  (11)  of  the  text. 


It  should  be  mentioned  that  the  final  results  may  readily  be  extended 
to  include  any  arbitrary  temporal  weighting  function  w(t).  This  is  most  easily 

(t-tb 

accomplished  by  replacing  the  uniform  weighting  function  n T in  (A-l)  by 
(t-t') 

w ' = *7(t-t')  II 1 c 1 > . This  substitution  carries  through  to  (A-10),  where 

A (r/ T)  must  be  replaced  by  l/T  [*?(r)  n ( r/T)  * y(r)  II  ( t/T)]  = 1/T 
[w(r ) ^w(r)],  Similarly,  \(r /T)  sine  2f  (T-  |r|)  in  (A-15)  must  be  replaced 

by  1/T[w(r)  e ^wr^<w(r)  e ^r]. 


Furthermore,  the  assumption  in  (A-16)  that  the  contribution  from 
integrals  involving  ^(x,  r)  is  negligible  in  comparison  to  the  contribution 
from  corresponding  integrals  involving  6(x,  r)  remains  valid  provided  that  f 
is  greater  than  the  width  of  the  spectral  window.  Thus,  in  most  practical 
situations,  the  only  modification  required  to  expressions  (A-17)  through  (A-19) 
is  that  A(  t/T)  in  (A-10)  be  replaced  by  1/T[w(r)  ^ w(r)J  . 


17/18 
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Appendix  B 


SIMPLIFICATION  OF  THE  BIAS  AND 
RANDOM  BEARING  ERROR  EXPRESSIONS 


In  this  appendix  we  wish  to  simplify  expressions  (A- 17)  through  (A- 19) 
or  (10)  and  (11)  in  the  text  by  using  the  following  two  additional  assumptions  of 
the  text: 

3.  The  signal  wavefront  is  a plane  wave  and  is  incident  on  the  array  at  an 
angle  60with  respect  to  the  normal  of  the  x-axis. 


4.  The  noise  is  spatially  uncorrelated  and  uniform. 


These  two  assumptions  allow  us  to  express  the  crosscorrelation 
between  any  two  points  in  the  space-time  aperture  as 


Rh(x,  T)  = Rg  (r-  e0  x)  +6(x)Rn(T), 


(B— 1) 


where  Rg  and  R are  the  signal  and  noise  autocorrelation  functions,  respec- 


tively; <5  (x)  is  the  familiar  impulse  function;  and£0is  defined  as  ~ ~°.  This 

effectively  uncouples  each  of  the  double  integrals  in  (A-17  ) through  (A-19) 
and  allows  each  to  be  expressed  as  the  sum  of  the  product  of  two  single 
integrals-one  integration  being  over  space  and  the  other  time.  To  illustrate 
this,  consider  an  integral  of  the  form 

I -JJ dxdr  4>(x,  T)  g(x),  (B-2) 

where  g(x)  depends  on  the  sensor  density  p(x).  Substituting  (A-10)  and  (B-l) 
into  (B-2)  yields 


I=//dxdrR  (r-  £0  x)  A (r/T)  e'iaT  ei^x  g(x) 


-iur  iu>£x 


Jf  dxd  T 6(x)  Rn  ( r)  A(  r/T ) e_1  ^ e 


(B-3) 


Letting  T_  - T-  £0x  in  the  first  double  integral,  we  obtain 


rr  , (V  ^°X)  -Mr  +{oX)  iwf. 

1=//  dxdriRg(rl)  A e 1 e *' 


(B-4) 


//dxdr  Rfl  (r)  A (r/T)  e i""r,s(x)  g(x)  eiw^x  . 
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Now,  the  maximum  value  that  x can  attain  is  ± L/2,  where  L is  the  total 
length  of  the  array.  Thus,  provided  T >>  L/c, 


A 


(Tj  + i0x) 
T 


« A(T1/T)  • 


Using 

/ dxft  (x)  elu"^X  g(x)  = g(0) 


reduces  (B-4) to 

I=/dTR  (y)  A(r/T)  e"iwT  •/  dx  g(x)  eiw(*_*°  )x 
' s * 

+ g(0)  «/dTRn(r)  A(T/T)  e~iujT. 


(B-5) 


An  equivalent  expression  for  the  above  integral  in  terms  of  the  signal  and 

noise  autospectral  densities,  G (f)  and  G (f),  can  be  obtained  by  noting  that 

s n 


1/T/drR  (r)  A(t/T)  e-luT  = f daG  (a)  sine2  T(d-f) 


= G (f)  * sine2  Tf  = S(f) 
s 


(B-6) 


and 


l/T/dTRn(T)  A( f/T)  e~iu,T=  J doGn(f)  sine2  T(o-f) 


= Gn(f)  * sine  Tf  = N(f)  , 


(B-7) 


where  S(f)  and  N(f)  are  the  respective  signal  and  noise  power  within  the 
processing  band  centered  at  frequency  f.  Using  these  definitions  and 
evaluating  the  integral  at  £ = £0  , we  obtain 


//dxd  r<i>(x,  r)  g(x)  = T S(f)/g(x)dx  + T N(f)  g(0).  (B- 

S <0 


8) 


Expression  (B-8)  may  now  be  used  to  evaluate  each  of  the  double 
integrals  in  (A-17)  through  (A-19).  Accordingly,  one  obtains  (in  shorthand 
notation) 


20 


where  ( ) denotes  that  the  resulting  function  is  to  be  evaluated  at  x = 0. 

The  fact  that  E(y.')  | ^ =0  implies  that  the  estimator  is  unbiased  under  the 

assumed  conditions.  0 


In  a lengthy  but  straightforward  operation,  expression  (B-9)  may  be 
expanded  and  rewritten  as 

E(y'2)  = u>2  {SN  [2/pdxy*x2  (p  * p2)  dx -/ p2  dx/x2  (p  * p)  dx] 


Notice  that  the  coefficients  of  S have  canceled 


Finally,  substitution  of  (B-10)  through  (B-12)  into  (6)  in  the  text 


yields 


£ (estimator  is  unbiased) 
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